Lyapunov exponents of heavy particles in turbulence 
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Lyapunov exponents of heavy particles and tracers advected by homogeneous and isotropic turbu- 
lent flows are investigated by means of direct numerical simulations. For large values of the Stokes 
number, the main effect of inertia is to reduce the chaoticity with respect to fluid tracers. Conversely, 
for small inertia, a counter-intuitive increase of the first Lyapunov exponent is observed. The flow 
intermittency is found to induce a Reynolds number dependency for the statistics of the finite time 
Lyapunov exponents of tracers. Such intermittency effects are found to persist at increasing inertia. 



Impurities suspended in incompressible flows are rele- 
vant to several physical processes, such as spray combus- 
tion, raindrop formation and transport of pollutantsiS*^ 
These particles, being typically of finite size and heav- 
ier than the ambient fluid, cannot be modeled as sim- 
ple tracers. They indeed posses inertia, which is respon- 
sible for the spontaneous generation of strong inhomo- 
geneities in their spatial distribution^i^ In a turbulent 
flow, their clustering is more efficient at small scales, be- 
low the dissipative scale, where the fluid velocity field 
is smooth. The intensity of particle clustering is there 
related to the statistics of the Lyapunov exponents of 
particle trajectoriesii The behavior of the Lyapunov ex- 
ponents of inertial particles and the relation with parti- 
cle clustering was recently investigated in random short- 
correlated flowSi^iSiiS 

In this Letter, be means of high-resolution direct nu- 
merical simulations, we investigate the Lyapunov spec- 
tra of inertial particles, varying their response time 
and the Reynolds number of the carrier turbulent flow. 
For Ts larger than the Kolmogorov-scale turnover time 
Tjj , the presence of inertia results in a generic reduction 
of chaoticity: the leading Lyapunov Exponent is smaller 
than the one of tracers and the strongest chaotic fluctu- 
ations become less probable. Remarkably, for < r^, 
an increase of chaoticity is observed; this effect can be 



understood in terms of the preferential concentration of 
particles in the high-strain regions of the flow. For fluid 
tracers, we observe intermittency effects on the statis- 
tics of chaotic fluctuations: the Reynolds-number depen- 
dence deviates from the dimensional prediction. Such 
deviations are found to persist for particles with inertia. 

A small spherical particle of radius a with a density 
p much larger than the density po of the surrounding 
incompressible fluid, evolves according to the dynamicsii 
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V , V ^ [v — u{x{t),t)] , 



(1) 



where u is the fluid velocity at the location x of the 
particle that moves with velocity v, and = 2a} pj (Oi'po) 
is the Stokes time [y is the kinematic viscosity of the 
fluid). Equation holds when the flow surrounding the 
particle is a Stokes flow; this requires a <^ rj, r] being 
the Kolmogorov scale of the flow. The Stokes number 
is defined as St = Ts/rr,, where is the eddy turnover 
time associated to rj. Particles are assumed to behave 
passively: we neglect their feedback on the flow, which is 
justified for very diluted suspensionsi^^ For comparison, 
we also study the motion of neutral particles that follow 
the dynamics x = u{x{t),t), which corresponds to the 
limit Ts in Eq. (jTJ. 
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The incompressible fluid velocity field u{x, t) evolves 
according to the Navier-Stokes equations 

dtu + u-Vu= -{l/po)Vp + iyAu + f , V-m = 0, (2) 

where the forcing provides an external energy input at 
a rate e = {f ■ u). These equations are integrated 
on d = 3 dimensional grid of size N = 128, 256, 512 
with periodic boundary conditions, by means of a fully 
de-aliased pseudo-spectral parallel code with 2"'^ order 
Adams-Bashforth time-stepping. Energy is injected by 
keeping constant the spectral content of the two smallest 
wavenumber shells. Viscosity is chosen to resolve well 
the small-scale velocity {kmaxV — 1-7). The Reynolds 
numbers based on Taylor's micro-scale is in the range 
Rx e [65 : 185] (see Ref.[l3for further details). 

Once the fluid flow has reached a statistically station- 
ary state, particles are homogeneously seeded with initial 
velocities equal to the fluid velocity at their locations. We 
followed 33 sets of 2000 particles with Stokes numbers 
ranging from to 2.2 for a time ~ 200 after relaxation 
of transients. 

In order to compute the Lyapunov spectrum we follow 
along each particle trajectory the time evolution of 2 x d 
infinitesimal displacements in the position- velocity phase 
space obtained by linearizing the particle dynamics 
The infinitesimal volume , defined by j linear indepen- 
dent tangent vectors, grows in time with an exponential 
rate ELi7i(T) = (1/T) In [yJ(r)/W(0)] , which de- 
fines the finite-time Lyapunov exponents ji(T) (FTLE), 
also called stretching rates. They asymptotically con- 
verge to the Lyapunov exponents Xi = limT^oo 7i {T) , 
which by definition are labeled in decreasing order Ai > 
• • • ^ X.2d- To compute these exponents numerically, 
we make use at each time lag r,, of a standard tech- 
nique based on the orthonormalization of the infinites- 
imal displacements by a Gram-Schmidt procedure (see, 
e.g., Ref.Ill. 

In Fig. n we show the behavior of the first three Lya- 
punov exponents as a function of St for the largest value 
of the Reynolds number. These three exponents rule the 
time evolution of infinitesimal elements in the physical 
space. For the range of Stokes numbers investigated here, 
we observe X4 ^ X5 ^ Xq ~ —1/ts, signaling the relax- 
ation of particle velocities to the fluid. 

The largest Lyapunov exponent Ai measures the 
chaotic separation of particle trajectories. To understand 
how chaoticity is affected by inertia, two mechanisms 
have to be considered. First, inertial particles have a 
delay on the fluid motion; this means that their velocity 
is approximately given by that of tracers with a time fil- 
tering over a time window of size r, . This effect weakens 
chaoticity. Second, heavy particles are ejected from per- 
sistent vortical structures and concentrate in high-strain 
regions. Since these portions of the flow are characterized 
by larger stretching rates, the chaoticity of particle tra- 
jectories is increased with respect to tracers that homoge- 
neously visit all the regions. As emphasized in the inset 
of Fig. n the latter effect dominates for St < 1 where Ai 
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FIG. 1: Lyapunov exponents Ai for i — 1 (triangles), i = 2 
(crosses) and i = 3 (circles) as a function of Stokes number. 
Rx ~ 185. In the inset we show the relative growth of the 
first Lyapunov exponent Xi{St)/ Xi{0) occurring at small St. 



is larger than the Lyapunov exponent of the fluid trac- 
ers {St — 0). This is not be predicted from analytical 
and numerical studies done in white-in-time random ve- 
locity fieldsiSiiS such flows clearly possess no persistent 
structures. Conversely, at sufficiently large St, the Lya- 
punov exponent decreases: preferential concentration is 
then negligible and the time-filtering approximation be- 
comes relevant. For such a large inertia, the white-in- 
time models apply and indeed predict a decrease of Ai 
as a function of St. Note that the competition between 
filtering and preferential concentration described above 
also enters in the distribution of particle acceleration>i^ 
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FIG. 2: Lyapunov dimension dx (squares) and volume growth 
rate A = X^i^i Xi (circles) as a function of Stokes number 
{R\ — 185). Inset: PDF of the finite-time volume growth 
rate r(r) for T ^ 80 r,,. 

As observed in Fig. ^ the time evolution of infinites- 
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imal surfaces is also affected by these two mechanisms. 
Indeed, at varying St, the second Lyapunov exponent A2 
displays a behavior qualitatively similar to that of Ai. 
For the tracers {St = 0), incompressibility of the flow 
implies Ai + A2 + A3 = 0. Since for a time reversible 
dynamics one has A2 ~ 0, the ratio A2/A1 is a measure 
of the irreversibility of the dynamics. Previous numer- 
ical investigations at moderate Reynolds numbcrsi^ in- 
dicate A2/A1 ~ 0.25; our simulations indicate A2/A1 ~ 
(0.28 ± 0.02). This irreversibility stems from the fact 
that the Navier-Stokes equation itself is not invariant 
with respect to time reversal. 
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FIG. 3: Cramer function S(7i) for fluid tracers for for R\ — 
65 (squares), 105 (circles) and 185 (triangles). Inset: Lya- 



punov exponent Ai (crosses) and reduced variance /i 
(diamonds) as a function of R\. 
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The behavior of A3 as a function of St markedly dif- 
fers from that of the first two exponents. Due to the 
dissipative nature of the inertial particle dynamics, vol- 
umes in physical space are not conserved. Indeed, the 
volume growth rate, defined as A = Ai -I- A2 -I- A3, which 
identically vanishes for fluid tracers, is negative for all 
Stokes numbers in the range < St < 1.72 (see Fig.|2l. 
This means that all volumes from physical space contract 
to zero at large times. Such clustering, which happens 
at scales much smaller than rj, is a consequence of the 
convergence of particle trajectories toward a dynamically 
evolving (multi)fractal setii The fractal dimension of this 
attractor can be estimated by means of the Kaplan- Yorkc 
(or Lyapunov) dimensioni^ as d\ = J + X^iLi K/\^.J+i\ 
where J is the largest integer such that X^iLi > 0. 
The fractal dimension of inertial particles is represented 
in Fig. 121 as a function of St. The minimum at St w 0.5 
corresponds to maximal clustering. For St > 1.72 the 
fractal dimension becomes greater than d ~ 3, indicat- 
ing that the spatial distribution of particles is not frac- 
tal anymore. For 5*^ ~ 1.72 the volume growth rate A 
vanishes, meaning that the dynamics of such particles 
preserve volumes on average. However, the finite-time 



volume growth rate F = 71 -f- 72 -f 73 experiences large 
fluctuations, as shown in the inset of Fig. |21 As a result, 
strong local inhomogeneities are present in the particle 
concentration also at large values of St. 

We now turn to the study of the Lyapunov exponent 
dependence on the Reynolds number of the flow. The 
inset of Fig. 13 shows the first Lyapunov exponent for 
neutral particles (i.e. St = 0) as a function of R\. Since 
Ai is a small-scale turbulent quantity with the dimension 
of an inverse time, one expects Ait^ ~ const and thus 
Ai oc R\. However, as a consequence of the intermittent 
fluctuations of the velocity gradients in turbulent flows, 
one can predict an anomalous dependence on Reynolds 
number Ai cx R" with a < 1—. This implies that Air,, 
decreases with Reynolds, as indeed confirmed by our sim- 
ulations. 




FIG. 4: Mean value Ai (small symbols) and reduced variance 
fi (large symbols) of the FTLE as a function of the Stokes 
number, for Rx = 65 (squares), Rx = 105 (circles) and R\ = 
185 (triangle). Inset: Cramer function 5(71) for Rx ~ 185, 
and various values of the Stokes number: St = (triangles), 
St = 0.32 (diamonds) and St = 2.19 (stars). 

Intermittency is actually expected to affect the whole 
probability distribution function (PDF) of the largest 
finite-time Lyapunov exponents 71 (T) . For T sufficiently 
large, the distribution of FTLE is expected to obey a 
large-deviation principle, i.e. pt{ii) oc exp(— TS'(7i)). 
The Cramer (or rate) function 5(71) is a non- negative 
concave function that vanishes at its minimum, attained 
for 71 = Ai. Small fluctuations occurring when I71 — 
All ^ T-1/2 are described by the central-limit theorem 
which amounts to approximating S{-) by a parabola in 
the vicinity of its minimum. The effect of intermittency 
on such small fluctuations can be measured from the vari- 
ance = ((71 — Ai)^) of the FTLE, or more particularly 
from the reduced variance ^ = Ta^ which measures the 
width of the Cramer function. As predicted in Ref.0, in- 
termittency is responsible for an anomalous dependency 
of ^ on the Reynolds number. More particularly ^r,, 
is expected to grow with R\. This tendency is qualita- 
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tively confirmed by our simulations as shown on the inset 
of Fig. 13 The signature of intermittency on the higher- 
order statistics of 71 can hardly be measured in a reliable 
way. Indeed, as shown in Fig. O the PDFs of 71 for the 
three Rx considered, once centered and normalized, al- 
most collapse for fluctuations as large as 3a. However 
it is still possible to observe a systematic deviation from 
the Gaussian distribution. Because of incompressibility, 
the left tail of the PDF is bounded by the constraint that 
7i > 0. It thus has to decrease faster than a Gaussian as 
indeed observed. The right part of the PDF is related to 
strong velocity gradients. Such events apparently lead to 
a tail which is fatter than Gaussian. Turbulent intermit- 
tency should therefore mainly affect the right tail. 

For inertial particles, intermittent corrections act in 
the same direction as for fluid tracers. Figure0|shows the 
behavior of Ai and of the reduced variance ^ as a func- 
tion of the Stokes number for various Reynolds numbers. 
As for tracers, for any given St, Air^ decreases while 
/iT,j increases with Rx- The inset of Fig. ^ shows, for 
Rx = 185, the Cramer function of the FTLE 71 (T) for 
both neutral particles {St = 0) and inertial particles with 
two different St. For St < 1, the whole distribution shifts 
to higher values, signaling the increased chaoticity. For 
St > 1, the distribution of FTLE shifts to lower values 
and fluctuations becomes less probable. The asymmetry 
observed in the PDF of tracer stretching rates decreases 
with inertia. At the largest values of St which are consid- 
ered, the Cramer function 71 becomes indistinguishable 



from a parabola. Finally, it is worth noticing that the 
dependence on the Reynolds number is less significant 
for the volume growth rate A, and thus for the fractal di- 
mension. Intermittency affects only weakly particle clus- 
tering. A clearer understanding of this issue requires to 
investigate larger values of the Reynolds number with 
longer statistics. This can hardly be achieved numeri- 
cally but could be done experimentally using standard 
box-counting techniques to estimate the fractal dimen- 
sions of the particle spatial distribution. 

We have seen that two mechanisms enter the dynam- 
ics of inertial particles: they concentrate in high-strain 
regions and they lag behind the fluid flow. Only in either 
the limit of small or large inertia, one of the two effect 
dominates the other. At present time, tackling analyti- 
cally the behavior of the largest Lyapunov exponent as a 
function of the Stokes number could only be done in these 
asymptoticS]^ For the range of Stokes numbers considered 
here, both effects compete and influence the Lyapunov 
exponents, preventing a complete analytical description. 
A numerical confirmation of the present theoretical pre- 
dictions would require larger computational resources. 
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